#Declarations
var N[y.end.sim, n.class], A[n.class,n.class,y.end.sim], rN[y.end.sim,n.class], f.neg[y.end.sim], f.pos[y.end.sim], f.con[y.end.sim], N.total[y.end.sim], phi[y.end.sim], r.sero.calf[y.end.sim], r.sero.yrlg[y.end.sim], r.sero.adult[y.end.sim]
#Data assingnment for initial conditions
data {
	mu.1 <-y.N[1]
	
	}
model{

#Priors 
#<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-
#Coeficients in recruitment regresion
#from analysis of YNP data



#Priors from Fuller paper
b[1] ~ dbeta(77,18.08)
b[2] ~ dbeta(36.5,20.54)
b[3] ~ dbeta(3.2, 11.36)


#Priors on pi and rho for serology test, from meta analysis
#alpha_pi ~ dgamma(10.7,7.8)
#beta_pi ~ dgamma(11.7,1.6)
#pi ~ dbeta(alpha_pi,beta_pi)
#
#alpha_rho ~ dgamma(7.8,7.8)
#beta_rho ~ dgamma(5.5,.19)
#rho ~ dbeta(alpha_rho,beta_rho)
pi ~ dbeta(476,42)
rho ~ dbeta(403,15467)


#shape parameters for gamma distribution of census standard deviations

alpha.var ~ dgamma(.001,.001)
nu.var ~ dgamma(.001,.001)

#prior on sighting probability
p.sight ~ dbeta(370,9.50)


#prior on vaccination efficacy


epsilon ~ dbeta(25+1, 74-25+1)

#Survival
#Juvenile
p[1] ~ dunif(0,1)
#Adult. Mean and sd of survival taken from YNP white binder
#mean<-.96, s[d <- .03, converted to shape parameters for beta
p[2] ~ dbeta(40, 1.667) 
p[3] ~ dunif(0,1)

#weakly informative prior on offspring sex ratio (mean <- .5, sd <- .05)
m ~ dbeta(49.5,49.5)

#disease parameters

beta ~ dunif(0,50)

v ~ dbeta(y.shapes.v[1],y.shapes.v[2])

psi ~ dbeta(y.shapes.psi[1], y.shapes.psi[2])

#probability of seropositive in removals
for(t in 1:y.end.sim){
	r.sero.calf[t] ~ dbeta(1,1)
	r.sero.yrlg[t] ~ dbeta(1,1)
	r.sero.adult[t] ~ dbeta(1,1)
	}

#stochasticity in process model


	sigma.p ~ dunif(0,5)
	tau.p <- 1/(sigma.p^2)	

#End of priors
#<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-

###inital conditions
#Initial conditions
for(j in 1:n.class){
	shape.d[j] <- 1
}
comp.1 ~ ddirch(shape.d[1:n.class])
tau.1 <- 1/var.N.obs
N.1 ~ dnorm(mu.1,tau.1)
for(j in 1:n.class){
	N.start[j]<- N.1*comp.1[j] 
}	

N[1,1:n.class] <- N.start[1:n.class]


#inital conditions for dervied quantities and latent states
# total population size

f.neg[1] <- b[1]
f.pos[1] <- b[2]
f.con[1]  <- b[3]

#z.calf[1] ~ dunif(0,100)  #These are needed to allow initial values for t > 1
#z.yrlg[1] ~ dunif(0,100)
#z.adult[1] ~ dunif(0,100)

diff.pos <- f.neg[1]-f.pos[1]
diff.con <- f.neg[1]-f.con[1]
N.total[1] <-sum(N[1,1:n.class])

#prevalance by age class
	pv.cow[1] <-  sum(N[1,6:7])/(sum(N[1,6:7])+N[1,3])
	pv.yrl[1] <- N[1,5] /(N[1,2] + N[1,5])
	pv.calf[1] <- N[1,4]/(N[1,4]+N[1,1]) 
#Population compostion 

	calf.ratio[1] <- (N[1,1] + N[1,4]) / N.total[1] 
	cow.ratio[1] <- (N[1,3] + sum(N[1,6:7]))/ N.total[1] 
	yrl.ratio[1] <- (N[1,2] + N[1,5])/ N.total[1]
	bull.ratio[1] <- N[1,8]/N.total[1]
	
#Probability of transmission

#Probability of transmission

z<-1
phi[1] <- 1-exp(-beta * (p[2]*N[1,6] + p[2]*psi*N[1,7])   / (1-m)*(p[1]*(N[1,1]+N[1,4]) + p[2]*(sum(N[1,2:3])+p[2]*sum(N[1,5:7]) ))^z)



#proportion of adult cows that are infectious
in.per[1] <- (N[1,6]/(sum(N[1,6:7]) + N[1,3]))
#propotion of sero-postive adult cows that are infectious
in.per.pos[1] <- (N[1,6])/(sum(N[1,6:7]))

#Force of infection
F[1] <- 0 #placeholder for first value in vector
#effective reproductive number
Re[1] <-0 #placeholder for first value in vector
#proportion of new infections by vertical transmission

v.b[1] <- p[2]* v *(f.con[1]*N[1,6] + f.pos[1]*N[1,7])
vert.per[1] <- v.b[1] / sum(N[1,4:6])


#estimate proportion of each age class in removals that is seropositive
for(i in 1:length(y.remove.index)){
	
	#z.calf[y.remove.index[i]] ~ dbin(r.sero.calf[y.remove.index[i]],y.rsero.calf[y.remove.index[i],2])
#	y.rsero.calf[y.remove.index[i],1] ~ dpois(pi*z.calf[y.remove.index[i]] + rho*(y.rsero.calf[y.remove.index[i],2] - z.calf[y.remove.index[i]]))
	
	z.calf[i] ~ dbin(r.sero.calf[y.remove.index[i]],y.rsero.calf[y.remove.index[i],2])
	y.rsero.calf[y.remove.index[i],1] ~ dbin((pi*z.calf[i] + rho*(y.rsero.calf[y.remove.index[i],2] - z.calf[i]))/y.rsero.calf[y.remove.index[i],2],y.rsero.calf[y.remove.index[i],2])
#		
	
	#z.yrlg[y.remove.index[i]] ~ dbin(r.sero.yrlg[y.remove.index[i]],y.rsero.yrlg[y.remove.index[i],2])
	#y.rsero.yrlg[y.remove.index[i],1] ~ dpois(pi*z.yrlg[t] + rho*(y.rsero.yrlg[t,2] - z.yrlg[t]))
	z.yrlg[i] ~ dbin(r.sero.yrlg[y.remove.index[i]],y.rsero.yrlg[y.remove.index[i],2])
 	y.rsero.yrlg[y.remove.index[i],1] ~ dbin((pi*z.yrlg[i] + rho*(y.rsero.yrlg[y.remove.index[i],2] - z.yrlg[i]))/y.rsero.yrlg[y.remove.index[i],2],y.rsero.yrlg[y.remove.index[i],2])
#	
	# z.adult[t] ~ dbin(r.sero.adult[t],y.rsero.adult[t,2])
	# y.rsero.adult[t,1] ~ dpois(pi*z.adult[t] + rho*(y.rsero.adult[t,2] - z.adult[t]))
	
	z.adult[i] ~ dbin(r.sero.adult[y.remove.index[i]],y.rsero.adult[y.remove.index[i],2])
	y.rsero.adult[y.remove.index[i],1] ~ dbin((pi*z.adult[i] + rho*(y.rsero.adult[y.remove.index[i],2] - z.adult[i]))/y.rsero.adult[y.remove.index[i],2],y.rsero.adult[y.remove.index[i],2])
}

#constant for portion of the year elapsed before peak of removals
	
q<-3/4


#Process Model<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-
#create a survival vector for use in calculating probability of transmission

for(t in 2:y.end.sim){

##cacluate recrutiment rates for sero positive and sero negative females
#Recruitment parameters
f.neg[t] <- b[1]
f.pos[t] <- b[2]
f.con[t]  <- b[3]


#removals
#####Matrix version of removals
#estimate age / sex composition of removals
	r.comp[t,1:4]~ddirch(y.rage[t,1:4])
	#r.sero.calf[t] ~ dbeta(y.rsero.calf[t,1], y.rsero.calf[t,2])
	#r.sero.yrlg[t] ~ dbeta(y.rsero.yrlg[t,1], y.rsero.yrlg[t,2])
	#r.sero.adult[t] ~ dbeta(y.rsero.adult[t,1], y.rsero.adult[t,2])
	
	#probabilities that a removed animal is seropositive	
		
	#seronegtives removed
	#calves
	rN[t,1] <-y.removals[t]*r.comp[t,1]*(1-r.sero.calf[t])
	#yearlings
	rN[t,2] <-y.removals[t]*r.comp[t,2]*(1-r.sero.yrlg[t])
	#adults
	rN[t,3] <-y.removals[t]*r.comp[t,3]*(1-r.sero.adult[t])
	#######sero positives removed
	#calves
	rN[t,4] <-y.removals[t]*r.comp[t,1]*r.sero.calf[t]
	#yearlings
	rN[t,5] <-y.removals[t]*r.comp[t,2]*r.sero.yrlg[t]
	#removals of seropositive adult cows must be partitioned among 6, and 7
	rN[t,6] <-y.removals[t]*r.comp[t,3]*r.sero.adult[t] * N[t-1,6]/sum(N[t-1,6:7]) 	
	rN[t,7] <-y.removals[t]*r.comp[t,3]*r.sero.adult[t] * N[t-1,7]/sum(N[t-1,6:7])  	
	
	
	
	#bulls removed
	rN[t,8] <- y.removals[t]*r.comp[t,4]
	
#suvival probability based on removals and natural mortality
	s[t,1] <- max(0, p[1]-p[1]^(1-q)*rN[t,1]/N[t-1,1])
	s[t,2] <- max(0, p[2]-p[2]^(1-q)*rN[t,2]/N[t-1,2])
	s[t,3] <- max(0, p[2]-p[2]^(1-q)*rN[t,3]/N[t-1,3])
	s[t,4] <- max(0, p[1]-p[1]^(1-q)*rN[t,4]/N[t-1,4])
	s[t,5] <- max(0, p[2]-p[2]^(1-q)*rN[t,5]/N[t-1,5])
	s[t,6] <- max(0, p[2]-p[2]^(1-q)*rN[t,6]/N[t-1,6])
	s[t,7] <- max(0, p[2]-p[2]^(1-q)*rN[t,7]/N[t-1,7])
	s[t,8] <- max(0, p[3]-p[3]^(1-q)*rN[t,8]/N[t-1,8])


##frequency depdendent transmission
numer_phi[t]<-s[t,6]*N[t-1,6] + s[t,7]*psi*N[t-1,7]
denom_phi[t] <- ((1-m)*(s[t,1]*N[t-1,1]+s[t,4]*N[t-1,4]) + sum(s[t,2:3]*N[t-1,2:3]) + sum(s[t,5:7]*N[t-1,5:7]))^z

phi[t] <- 1-exp(-beta * numer_phi[t] /denom_phi[t])


  






#####Transition Matrix
#row 1
A[1,3,t]<-s[t,3]*f.neg[t]*(1-phi[t])
A[1,6,t]<-s[t,6]*f.con[t]*(1-v)*(1-phi[t])
A[1,7,t]<-s[t,7]* (f.pos[t] *(1-psi)*(1-phi[t])  +  f.pos[t] * psi*(1-v)*(1-phi[t]))


#row 2
A[2,1,t]<-s[t,1]*(1-m)*(1-phi[t])
#row 3
A[3,2,t]<-s[t,2]*(1-phi[t])
A[3,3,t]<-s[t,3]*(1-phi[t])

#row4
A[4,3,t]<-s[t,3]*f.neg[t]*phi[t]
A[4,6,t]<-s[t,6]*f.con[t]*(v+phi[t]-v*phi[t])
A[4,7,t]<-s[t,7]*(f.pos[t]*psi*(v+phi[t]-v*phi[t]) + f.pos[t]*(1-psi)*phi[t])

#row5
A[5,1,t]<-s[t,1]*phi[t]*(1-m)
A[5,4,t]<-s[t,4]*(1-m)

#row6
A[6,2,t]<-s[t,2]*phi[t]
A[6,3,t]<-s[t,3]*phi[t]
A[6,5,t]<-s[t,5]
A[6,7,t] <-0

#row 7
A[7,6,t]<-s[t,6]
A[7,7,t]<-s[t,7]

#row 8
A[8,1,t]<-s[t,1]*m
A[8,4,t]<-s[t,4]*m
A[8,8,t]<-s[t,8]






#estimates of mean of posterior distribution
mu[t,1:n.class] <- A[,,t]%*%N[t-1,1:n.class]

#get logs of deterministic predictions
for(j in 1:8){
	log_mu[t,j] <- log(max(1,mu[t,j]))
	}

log_N[t,1:8] ~ dmnorm(log_mu[t,1:8],tau.p*y.I) 

#exponentiate log of population size
for(j in 1:8){
	N[t,j] <- min(exp(log_N[t,j]),10000)
	#N[t,j] <- exp(log_N[t,j])
	}
N.total[t] <- sum(N[t,1:n.class])

#Population compostion 

	calf.ratio[t] <- (N[t,1] + N[t,4]) / N.total[t] #calves
	cow.ratio[t] <- (N[t,3] + sum(N[t,6:7])) / N.total[t]
 #adult females
 	yrl.ratio[t] <- (N[t,2] + N[t,5]) / N.total[t]

	bull.ratio[t]<-N[t,8]/N.total[t]
	
	

#sero-prevalace by age
	pv.cow[t] <-  sum(N[t,6:7])/(sum(N[t,6:7])+N[t,3])
	pv.yrl[t] <- N[t,5] /(N[t,2] + N[t,5])
	pv.calf[t] <- N[t,4]/(N[t,4]+N[t,1]) 
	
#proportion of adult cows that are infectious
in.per[t] <- (N[t,6]/(sum(N[t,6:7]) + N[t,3]))
#propotion of sero-postive adult cows that are infectious
in.per.pos[t] <- (N[t,6])/(sum(N[t,6:7]))
#proportion of new infections by vertical transmission

v.b[t] <- s[t,2]* v *(f.con[t]*N[t-1,6] + f.pos[1]*N[t-1,7])
vert.per[t] <- v.b[t] / sum(N[t,4:6])

#force of infection
F[t] <- phi[t] * sum(s[t,1:3]*N[t,1:3]) / sum(N[t-1,1:3])

#Effective reproductive ratio
Re[t] <- phi[t] * sum(s[t,1:3]*N[t,1:3]) / (N[t-1,6]+psi*N[t-1,7])


} #end process model


#averages across years
xcalf.ratio <- mean(calf.ratio[6:41])
xcow.ratio <- mean(cow.ratio[6:41])
xyrl.ratio <- mean(yrl.ratio[6:41])
xbull.ratio <- mean(bull.ratio[6:41])

xpv.cow <- mean(pv.cow[6:41])
xpv.yrl <- mean(pv.yrl[6:41])
xpv.calf <- mean(pv.calf[6:41])

xin.per <- mean(in.per[6:41])
xin.per.pos <- mean(in.per.pos[6:41])
xvert.per <- mean(vert.per[6:41])

xphi <- mean(phi[6:41])

#Data models<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-

#Likelihood for annual counts################################
##Estimate of observation error in census data

for(i in 1:length(y.N.varindex[])){
	#estimate variance based on observed variance
	y.N.var[y.N.varindex[i]] ~ dgamma(alpha.var,nu.var)
}

#Estimate of the mean observation variance

var.N.obs <- (alpha.var/nu.var)

for(t in 1:length(y.N.index)){
	#method of moments to get gamma shape parameters
	#from model estimates
	alpha[t]<- N.total[y.N.index[t]]^2/var.N.obs
	nu[t] <- N.total[t]/var.N.obs
	lambda[t] ~ dgamma(alpha[t],nu[t]) 
	lambda.cor[t]<-p.sight*lambda[t]
	y.N[t] ~ dpois(lambda.cor[t])
	
	#simulated data for posterior predictive check
	y.N.new[t] ~ dpois(p.sight*lambda[t])
	y.N.sq[t] <- (p.sight*lambda[t] - y.N[t])^2
	y.N.new.sq[t] <- (p.sight*lambda[t]-y.N.new[t])^2
}
	#Bayesian P value
	P.N.mean <- step(mean(y.N.new[]) - mean(y.N[]))
	P.N.sd <- step(sd(y.N.new[]) - sd(y.N[]))
	P.N.sq <- step(sum(y.N.new.sq[]) - sum( y.N.sq[]))
	
##  Likelihoods for population composition data ##############################################


#proportion of calves in population from aerial counts
for(t in 1:length(y.iratio.calf)){    #
	mu.calf[t] <- (N[y.iratio.calf[t],1] + N[y.iratio.calf[t],4]) / N.total[y.iratio.calf[t]]
#calcuate shape parameters for beta from data and estimates from external, hierachical model of ratio.calf data, y.sigma.calf[t] and ratio.calf[t] which are read in as data
	shape1.calf[t] <- max(.0001,(mu.calf[t]^2-mu.calf[t]^3-mu.calf[t]*y.sigma.calf[t]^2)/y.sigma.calf[t]^2)
	shape2.calf[t] <- max(.0001, (mu.calf[t]-2*mu.calf[t]^2+mu.calf[t]^3-y.sigma.calf[t]^2+mu.calf[t]*y.sigma.calf[t]^2)/y.sigma.calf[t]^2)
	y.ratio.calf[t] ~ dbeta(shape1.calf[t], shape2.calf[t])
	#simulate new data for posterior predictive checks
	y.ratio.calf.new[t] ~ dbeta(shape1.calf[t], shape2.calf[t])
	y.ratio.calf.sq[t] <- (mu.calf[t] - y.ratio.calf[t])^2
	y.ratio.calf.new.sq[t] <- (mu.calf[t] - y.ratio.calf.new[t])^2

} 
#Bayesian P values

P.ratio.calf.mean <- step(mean(y.ratio.calf.new[])-mean(y.ratio.calf[]))
	P.ratio.calf.sd <- step(sd(y.ratio.calf.new[])-sd(y.ratio.calf[]))
	P.ratio.calf.sq <- step(sum(y.ratio.calf.new.sq[]) - sum(y.ratio.calf.sq[]))


#proportions of yearling cows, adult cows, and non-yearling bulls in the population from ground counts
#for the first two years, the yearlings were lumped with the adult females. 
for(t in 1:2){
	alpha.p[t,1] <- bull.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,2] <- (cow.ratio[y.ground.index[t]] + yrl.ratio[y.ground.index[t]])*y.alpha_0[t]
	alpha.p[t,3] <- 0.0001*y.alpha_0[t]
	alpha.p[t,4] <- calf.ratio[y.ground.index[t]]*y.alpha_0[t]
	y.p.mu[t,1:4] ~ ddirch(alpha.p[t,1:4])
	}
for(t in 3:length(y.ground.index[])){
	alpha.p[t,1] <- bull.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,2] <- cow.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,3] <- yrl.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,4] <- calf.ratio[y.ground.index[t]]*y.alpha_0[t]
	y.p.mu[t,1:4] ~ ddirch(alpha.p[t,1:4])
	}






#Likelihoods for serology data####################################

#serology of calves:
for(j in 1:length(y.isero.calf[])){
		
	y.pos.calf[j] ~ dbin((N[y.isero.calf[j],4]*pi + N[y.isero.calf[j],1]*rho)/ (N[y.isero.calf[j],4]+N[y.isero.calf[j],1]), y.n.calf[j])
	
	#simulate data for posterior predictive check
	y.pos.calf.new[j] ~ dbin((N[y.isero.calf[j],4]*pi + N[y.isero.calf[j],1]*rho)/ (N[y.isero.calf[j],4]+N[y.isero.calf[j],1]), y.n.calf[j])

	y.pos.calf.sq[j] <- (y.pos.calf[j]-pv.calf[j])^2
	y.pos.calf.sq.new[j] <- (y.pos.calf.new[j] - pv.calf[j])^2
			
	}
	#Bayesian P value
	P.sero.mean[1] <- step(mean(y.pos.calf.new[])-mean(y.pos.calf[]))
	P.sero.sd[1] <- step(sd(y.pos.calf.new[])-sd(y.pos.calf[]))
	P.sero.sq[1] <- step(sum(y.pos.calf.sq.new[])-sum(y.pos.calf.sq[]))

	
#serology of yearlings
for(j in 1:length(y.isero.yrl[])){
	#sero-prevalance in calves
	
	y.pos.yrl[j] ~ dbin((N[y.isero.yrl[j],5]*pi + N[y.isero.yrl[j],2]*rho)/ (N[y.isero.yrl[j],5]+N[y.isero.yrl[j],2]), y.n.yrl[j])
	
	#simulate data for posterior predictive check
	y.pos.yrl.new[j]~ dbin((N[y.isero.yrl[j],5]*pi + N[y.isero.yrl[j],2]*rho)/ (N[y.isero.yrl[j],5]+N[y.isero.yrl[j],2]), y.n.yrl[j])

	y.pos.yrl.sq[j] <- (y.pos.yrl[j]-pv.yrl[j])^2
	y.pos.yrl.sq.new[j] <- (y.pos.yrl.new[j] - pv.yrl[j])^2
	}
	
	#Bayesian P value
	P.sero.mean[2] <- step(mean(y.pos.yrl.new[])-mean(y.pos.yrl[]))
	P.sero.sd[2] <- step(sd(y.pos.yrl.new[])-sd(y.pos.yrl[]))
	P.sero.sq[2] <- step(sum(y.pos.yrl.sq.new[])-sum(y.pos.yrl.sq[]))

	
#serology of adult cows
for(j in 1:length(y.isero.cow[])){
	y.pos.cow[j] ~ dbin((sum(N[y.isero.cow[j],6:7])*pi + N[y.isero.cow[j],3]*rho)/ (sum(N[y.isero.cow[j],6:7])+N[y.isero.cow[j],3]), y.n.cow[j])

	
	#Simulate data for posterior predictive check
	y.pos.cow.new[j] ~ dbin((sum(N[y.isero.cow[j],6:7])*pi + N[y.isero.cow[j],3]*rho)/ (sum(N[y.isero.cow[j],6:7])+N[y.isero.cow[j],3]), y.n.cow[j])
	y.pos.cow.sq[j] <- (y.pos.cow[j]-pv.cow[j])^2
	y.pos.cow.sq.new[j] <- (y.pos.cow.new[j] - pv.cow[j])^2

	}
	#Bayesian P value
	P.sero.mean[3] <- step(mean(y.pos.cow.new[])-mean(y.pos.cow[]))
	P.sero.sd[3] <- step(sd(y.pos.cow.new[])-sd(y.pos.cow[]))
	P.sero.sq[3] <- step(sum(y.pos.cow.sq.new[])-sum(y.pos.cow.sq[]))





##Forecasts
## Null model: do nothing


for(j in 1){
	n.s[1,1:8,j]<-N[y.end.sim,1:8]  #start with conditions at end of data
	n.s[1,9,j]<-0 
	n.s.total[1,j] <- N.total[y.end.sim]
	Re.s[1,j] <- Re[y.end.sim]
	F.s[1,j] <- F[y.end.sim]
	
	pv.cow.s[1,j] <-  pv.cow[y.end.sim]
	phi.s[1,j] <- phi[y.end.sim]
	calf.ratio.s[1,j] <- calf.ratio[y.end.sim] #calves
	cow.ratio.s[1,j] <- cow.ratio[y.end.sim]
	bull.ratio.s[1,j] <- bull.ratio[y.end.sim]
	
	delta.phi.s[1,j] <- 0
	#set placeholders for removals to allow output by setting trace monitor
	#set placeholders for removals to allow output by setting trace monitor
	r.ns[1,1,j]<-0
	r.ns[1,2,j]<-0  
	r.ns[1,3,j]<-0
	r.ns[1,4,j]<-0
	r.ns[1,5,j]<-0		
	r.ns[1,6,j]<-0
 	r.ns[1,7,j]<-0
	r.ns[1,8,j]<-0
	r.ns[1,9,j]<-0 #class for vacination in later runs

for(t in 2:y.scen.yrs+1){ #j indexes the scenario number
	
	#removal vector
	r.ns[t,1,j]<-0
	r.ns[t,2,j]<-0  
	r.ns[t,3,j]<-0
	r.ns[t,4,j]<-0
	#seropositives removed
	r.ns[t,5,j]<-0
	r.ns[t,6,j]<-0
	r.ns[t,7,j]<-0
	r.ns[t,8,j]<-0
	r.ns[t,9,j]<-0 #class for vacination in later runs
	
	   	
	#suvival probability based on removals and natural mortality
	S[t,1,j] <- p[1]
	S[t,2,j] <- p[2]
	S[t,3,j] <- p[2]
	S[t,4,j] <- p[1]
	S[t,5,j] <- p[2]
	S[t,6,j] <- p[2]
	S[t,7,j] <- p[2]
	S[t,8,j] <- p[3]
	S[t,9,j] <- p[2] #class for vacination in later runs


##frequency depdendent transmission
s_numer_phi[t,j]<-S[t,6,j]*n.s[t-1,6,j] + S[t,7,j]*psi*n.s[t-1,7,j]
s_denom_phi[t,j] <- ((1-m)*(S[t,1,j]*n.s[t-1,1,j]+S[t,4,j]*n.s[t-1,4,j]) + sum(S[t,2:3,j]*n.s[t-1,2:3,j]) + sum(S[t,5:7,j]*n.s[t-1,5:7,j]))^z
phi.s[t,j] <- 1-exp(-beta * s_numer_phi[t,j] /s_denom_phi[t,j])

#####Transition Matrix
#row 1
B[1,3,t,j]<-S[t,3,j]*f.neg[t]*(1-phi.s[t,j])
B[1,6,t,j]<-S[t,6,j]*f.con[t]*(1-v)*(1-phi.s[t,j])
B[1,7,t,j]<-S[t,7,j]* (f.pos[t] *(1-psi)*(1-phi.s[t,j])  +  f.pos[t] * psi*(1-v)*(1-phi.s[t,j]))
B[1,9,t,j]<-0

#row 2
B[2,1,t,j]<-S[t,1,j]*(1-m)*(1-phi.s[t,j])
#row 3
B[3,2,t,j]<-S[t,2,j]*(1-phi.s[t,j])
B[3,3,t,j]<-S[t,3,j]*(1-phi.s[t,j])

#row4
B[4,3,t,j]<-S[t,3,j]*f.neg[t]*phi.s[t,j]
B[4,6,t,j]<-S[t,6,j]*f.con[t]*(v+phi.s[t,j]-v*phi.s[t,j])
B[4,7,t,j]<-S[t,7,j]*(f.pos[t]*psi*(v+phi.s[t,j]-v*phi.s[t,j]) + f.pos[t]*(1-psi)*phi.s[t,j])
B[4,9,t,j]<-0
#row5
B[5,1,t,j]<-S[t,1,j]*phi.s[t,j]*(1-m)
B[5,4,t,j]<-S[t,4,j]*(1-m)

#row6
B[6,2,t,j]<-S[t,2,j]*phi.s[t,j]
B[6,3,t,j]<-S[t,3,j]*phi.s[t,j]
B[6,5,t,j]<-S[t,5,j]
B[6,7,t,j] <-0

#row 7
B[7,6,t,j]<-S[t,6,j]
B[7,7,t,j]<-S[t,7,j]

#row 8
B[8,1,t,j]<-S[t,1,j]*m
B[8,4,t,j]<-S[t,4,j]*m
B[8,8,t,j]<-S[t,8,j]


#row 9
B[9,9,t,j]<-S[t,9,j]
	
#estimates of mean of posterior distribution
mu.s[t,1:9,j] <- B[1:9,1:9,t,j]%*%n.s[t-1,1:9,j]

#get logs of deterministic predictions
for(k in 1:9){
	log_mu.s[t,k,j] <- log(max(1,mu.s[t,k,j]))
	}

log_n.s[t,1:9,j] ~ dmnorm(log_mu.s[t,1:9,j],tau.p*y.I2) 

#exponentiate log of population size
for(k in 1:9){
	n.s[t,k,j] <- min(exp(log_n.s[t,k,j]),10000)
	}

	
#derived quantities
Re.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / (n.s[t-1,6,j]+psi*n.s[t-1,7,j])
F.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / sum(n.s[t-1,1:3,j]) 
pv.cow.s[t,j] <-  sum(n.s[t,6:7,j])/(sum(n.s[t,6:7,j])+n.s[t,3,j])
n.s.total[t,j] <- sum(n.s[t,1:8,j])
cow.ratio.s[t,j] <- (n.s[t,3,j] + sum(n.s[t,6:7,j])) / n.s.total[t,j] #cows
calf.ratio.s[t,j] <- (n.s[t,1,j] + n.s[t,4,j]) / n.s.total[t,j] #calves
bull.ratio.s[t,j] <- n.s[t,8,j] / n.s.total[t,j] #bulls

delta.phi.s[t,j] <-phi.s[t,j] - phi[y.end.sim]
r.ns.total[t,j] <- sum(r.ns[t,1:8,j])

		} # end of t

x.r.ns[j] <- mean(r.ns.total[2:6,j])		


		
		} #end of  j

####remove seropositives##########################
##
for(j in 2){
	n.s[1,1:8,j]<-N[y.end.sim,1:8]  #start with conditions at end of data
	n.s[1,9,j]<-0 
	n.s.total[1,j] <- N.total[y.end.sim]
	Re.s[1,j] <- Re[y.end.sim]
	F.s[1,j] <- F[y.end.sim]
	
	pv.cow.s[1,j] <-  pv.cow[y.end.sim]
	phi.s[1,j] <- phi[y.end.sim]
	calf.ratio.s[1,j] <- calf.ratio[y.end.sim] #calves
	cow.ratio.s[1,j] <- cow.ratio[y.end.sim]
	bull.ratio.s[1,j] <- bull.ratio[y.end.sim]
	delta.phi.s[1,j] <- 0
	#set placeholders for removals to allow output by setting trace monitor
	r.ns[1,1,j]<-0
	r.ns[1,2,j]<-0  
	r.ns[1,3,j]<-0
	r.ns[1,4,j]<-0
	r.ns[1,5,j]<-0		
	r.ns[1,6,j]<-0
 	r.ns[1,7,j]<-0
	r.ns[1,8,j]<-0
	r.ns[1,9,j]<-0 #class for vacination in later runs
for(t in 2:y.scen.yrs+1){ #j indexes the scenario number
	
###This sets removals for all secnarios--producing the same removal number for each time
	prop.remove[t] ~ dbeta(y.a,y.b)
	
	r.remove[t,j] <-min(prop.remove[t]*n.s.total[t-1,j],1500)

	

	#removal vector
	r.ns[t,1,j]<-0
	r.ns[t,2,j]<-0  
	r.ns[t,3,j]<-0
	r.ns[t,4,j]<-0
	#seropositives removed--indexed to j = 1 to assure same number of removals
	r.ns[t,5,j]<-r.remove[t,j]*n.s[t-1,5,j]/n.s.total[t-1,j] 
	r.ns[t,6,j]<-r.remove[t,j]*n.s[t-1,6,j]/n.s.total[t-1,j]
	r.ns[t,7,j]<-r.remove[t,j]*n.s[t-1,7,j]/n.s.total[t-1,j]	
	r.ns[t,8,j]<-0
	r.ns[t,9,j]<-0 #class for vacination in later runs
	
	
	
	#suvival probability based on removals and natural mortality
	S[t,1,j] <- p[1]
	S[t,2,j] <- p[2]
	S[t,3,j] <- p[2]
	S[t,4,j] <- p[1]
	S[t,5,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,5,j]/n.s[t-1,5,j])
	S[t,6,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,6,j]/n.s[t-1,6,j])
	S[t,7,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,7,j]/n.s[t-1,7,j])
	S[t,8,j] <- p[3]
	S[t,9,j] <- p[3]
	
		

#
##frequency depdendent transmission
s_numer_phi[t,j]<-S[t,6,j]*n.s[t-1,6,j] + S[t,7,j]*psi*n.s[t-1,7,j]
s_denom_phi[t,j] <- ((1-m)*(S[t,1,j]*n.s[t-1,1,j]+S[t,4,j]*n.s[t-1,4,j]) + sum(S[t,2:3,j]*n.s[t-1,2:3,j]) + sum(S[t,5:7,j]*n.s[t-1,5:7,j]))^z
phi.s[t,j] <- 1-exp(-beta * s_numer_phi[t,j] /s_denom_phi[t,j])




#####Transition Matrix
#row 1
B[1,3,t,j]<-S[t,3,j]*f.neg[t]*(1-phi.s[t,j])
B[1,6,t,j]<-S[t,6,j]*f.con[t]*(1-v)*(1-phi.s[t,j])
B[1,7,t,j]<-S[t,7,j]* (f.pos[t] *(1-psi)*(1-phi.s[t,j])  +  f.pos[t] * psi*(1-v)*(1-phi.s[t,j]))
B[1,9,t,j]<-0

#row 2
B[2,1,t,j]<-S[t,1,j]*(1-m)*(1-phi.s[t,j])
#row 3
B[3,2,t,j]<-S[t,2,j]*(1-phi.s[t,j])
B[3,3,t,j]<-S[t,3,j]*(1-phi.s[t,j])

#row4
B[4,3,t,j]<-S[t,3,j]*f.neg[t]*phi.s[t,j]
B[4,6,t,j]<-S[t,6,j]*f.con[t]*(v+phi.s[t,j]-v*phi.s[t,j])
B[4,7,t,j]<-S[t,7,j]*(f.pos[t]*psi*(v+phi.s[t,j]-v*phi.s[t,j]) + f.pos[t]*(1-psi)*phi.s[t,j])
B[4,9,t,j]<-0
#row5
B[5,1,t,j]<-S[t,1,j]*phi.s[t,j]*(1-m)
B[5,4,t,j]<-S[t,4,j]*(1-m)

#row6
B[6,2,t,j]<-S[t,2,j]*phi.s[t,j]
B[6,3,t,j]<-S[t,3,j]*phi.s[t,j]
B[6,5,t,j]<-S[t,5,j]
B[6,7,t,j] <-0

#row 7
B[7,6,t,j]<-S[t,6,j]
B[7,7,t,j]<-S[t,7,j]

#row 8
B[8,1,t,j]<-S[t,1,j]*m
B[8,4,t,j]<-S[t,4,j]*m
B[8,8,t,j]<-S[t,8,j]


#row 9
B[9,9,t,j]<-S[t,9,j]
	
#estimates of mean of posterior distribution
mu.s[t,1:9,j] <- B[1:9,1:9,t,j]%*%n.s[t-1,1:9,j]

#get logs of deterministic predictions
for(k in 1:9){
	log_mu.s[t,k,j] <- log(max(1,mu.s[t,k,j]))
	}

log_n.s[t,1:9,j] ~ dmnorm(log_mu.s[t,1:9,j],tau.p*y.I2) 

#exponentiate log of population size
for(k in 1:9){
	n.s[t,k,j] <- min(exp(log_n.s[t,k,j]),10000)
	}

	
#derived quantities
Re.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / (n.s[t-1,6,j]+psi*n.s[t-1,7,j])
F.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / sum(n.s[t-1,1:3,j]) 
pv.cow.s[t,j] <-  sum(n.s[t,6:7,j])/(sum(n.s[t,6:7,j])+n.s[t,3,j])
n.s.total[t,j] <- sum(n.s[t,1:8,j])
cow.ratio.s[t,j] <- (n.s[t,3,j] + sum(n.s[t,6:7,j])) / n.s.total[t,j] #cows
calf.ratio.s[t,j] <- (n.s[t,1,j] + n.s[t,4,j]) / n.s.total[t,j] #calves
bull.ratio.s[t,j] <- n.s[t,8,j] / n.s.total[t,j] #bulls

delta.phi.s[t,j] <-phi.s[t,j] - phi[y.end.sim]

r.ns.total[t,j] <- sum(r.ns[t,1:8,j])

		} #end of t
		
x.r.ns[j] <- mean(r.ns.total[2:6,j])	
				
		} #end of j
###remove seronegatives##########################
###
for(j in 3){
	n.s[1,1:8,j]<-N[y.end.sim,1:8]  #start with conditions at end of data
	n.s[1,9,j]<-0 
	n.s.total[1,j] <- N.total[y.end.sim]
	Re.s[1,j] <- Re[y.end.sim]
	F.s[1,j] <- F[y.end.sim]
	
	pv.cow.s[1,j] <-  pv.cow[y.end.sim]
	phi.s[1,j] <- phi[y.end.sim]
	calf.ratio.s[1,j] <- calf.ratio[y.end.sim] #calves
	cow.ratio.s[1,j] <- cow.ratio[y.end.sim]
	bull.ratio.s[1,j] <- bull.ratio[y.end.sim]
	delta.phi.s[1,j] <- 0
	r.ns[1,1,j]<-0
	r.ns[1,2,j]<-0  
	r.ns[1,3,j]<-0
	r.ns[1,4,j]<-0
	r.ns[1,5,j]<-0		
	r.ns[1,6,j]<-0
 	r.ns[1,7,j]<-0
	r.ns[1,8,j]<-0
	r.ns[1,9,j]<-0
	
for(t in 2:y.scen.yrs+1){ #j indexes the scenario number
	
	#get random removal total, also used with remove seronegative and vaccination scenarios
	#a.mu[t,j] <-y.remove.mean
#	a.var[t,j] <-y.remove.var
#	l[t,j] ~ dgamma(a.mu[t,j]^2/a.var[t,j],a.mu[t,j]/a.var[t,j])
#	r.remove[t,j] ~ dpois(l[t,j])



	r.remove[t,j] <-min(prop.remove[t]*n.s.total[t-1,j],1500)
	
	#removal vector
	r.ns[t,1,j]<-0
	r.ns[t,2,j]<-r.remove[t,j]*n.s[t-1,2,j]/n.s.total[t-1,j] 	
	r.ns[t,3,j]<-r.remove[t,j]*n.s[t-1,3,j]/n.s.total[t-1,j] 
	r.ns[t,4,j]<-0
	#seropositives removed
	r.ns[t,5,j]<-0
	r.ns[t,6,j]<-0
	r.ns[t,7,j]<-0
	r.ns[t,8,j]<-0
	r.ns[t,9,j]<-0 #class for vacination in later runs
	
	
	S[t,1,j] <- p[1]
	S[t,2,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,2,j]/n.s[t-1,2,j])
	S[t,3,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,3,j]/n.s[t-1,3,j])
	S[t,4,j] <- p[1]
	S[t,5,j] <- p[2]
	S[t,6,j] <- p[2]	
	S[t,7,j] <- p[2]
	S[t,8,j] <- p[3]
	S[t,9,j] <- p[3]
	
	

#

#
##frequency depdendent transmission
s_numer_phi[t,j]<-S[t,6,j]*n.s[t-1,6,j] + S[t,7,j]*psi*n.s[t-1,7,j]
s_denom_phi[t,j] <- ((1-m)*(S[t,1,j]*n.s[t-1,1,j]+S[t,4,j]*n.s[t-1,4,j]) + sum(S[t,2:3,j]*n.s[t-1,2:3,j]) + sum(S[t,5:7,j]*n.s[t-1,5:7,j]))^z
phi.s[t,j] <- 1-exp(-beta * s_numer_phi[t,j] /s_denom_phi[t,j])
	

#####Transition Matrix
#row 1
B[1,3,t,j]<-S[t,3,j]*f.neg[t]*(1-phi.s[t,j])
B[1,6,t,j]<-S[t,6,j]*f.con[t]*(1-v)*(1-phi.s[t,j])
B[1,7,t,j]<-S[t,7,j]* (f.pos[t] *(1-psi)*(1-phi.s[t,j])  +  f.pos[t] * psi*(1-v)*(1-phi.s[t,j]))
B[1,9,t,j]<-0

#row 2
B[2,1,t,j]<-S[t,1,j]*(1-m)*(1-phi.s[t,j])
#row 3
B[3,2,t,j]<-S[t,2,j]*(1-phi.s[t,j])
B[3,3,t,j]<-S[t,3,j]*(1-phi.s[t,j])

#row4
B[4,3,t,j]<-S[t,3,j]*f.neg[t]*phi.s[t,j]
B[4,6,t,j]<-S[t,6,j]*f.con[t]*(v+phi.s[t,j]-v*phi.s[t,j])
B[4,7,t,j]<-S[t,7,j]*(f.pos[t]*psi*(v+phi.s[t,j]-v*phi.s[t,j]) + f.pos[t]*(1-psi)*phi.s[t,j])
B[4,9,t,j]<-0
#row5
B[5,1,t,j]<-S[t,1,j]*phi.s[t,j]*(1-m)
B[5,4,t,j]<-S[t,4,j]*(1-m)

#row6
B[6,2,t,j]<-S[t,2,j]*phi.s[t,j]
B[6,3,t,j]<-S[t,3,j]*phi.s[t,j]
B[6,5,t,j]<-S[t,5,j]
B[6,7,t,j] <-0

#row 7
B[7,6,t,j]<-S[t,6,j]
B[7,7,t,j]<-S[t,7,j]

#row 8
B[8,1,t,j]<-S[t,1,j]*m
B[8,4,t,j]<-S[t,4,j]*m
B[8,8,t,j]<-S[t,8,j]


#row 9
B[9,9,t,j]<-S[t,9,j]
	
#estimates of mean of posterior distribution
mu.s[t,1:9,j] <- B[1:9,1:9,t,j]%*%n.s[t-1,1:9,j]

#get logs of deterministic predictions
for(k in 1:9){
	log_mu.s[t,k,j] <- log(max(1,mu.s[t,k,j]))
	}

log_n.s[t,1:9,j] ~ dmnorm(log_mu.s[t,1:9,j],tau.p*y.I2) 

#exponentiate log of population size
for(k in 1:9){
	n.s[t,k,j] <- min(exp(log_n.s[t,k,j]),10000)
	}

	
#derived quantities
Re.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / (n.s[t-1,6,j]+psi*n.s[t-1,7,j])
F.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / sum(n.s[t-1,1:3,j]) 
pv.cow.s[t,j] <-  sum(n.s[t,6:7,j])/(sum(n.s[t,6:7,j])+n.s[t,3,j])
n.s.total[t,j] <- sum(n.s[t,1:8,j])
cow.ratio.s[t,j] <- (n.s[t,3,j] + sum(n.s[t,6:7,j])) / n.s.total[t,j] #cows
calf.ratio.s[t,j] <- (n.s[t,1,j] + n.s[t,4,j]) / n.s.total[t,j] #calves
bull.ratio.s[t,j] <- n.s[t,8,j] / n.s.total[t,j] #bulls

delta.phi.s[t,j] <-phi.s[t,j] - phi[y.end.sim]

r.ns.total[t,j] <- sum(r.ns[t,1:8,j])

		} #end of t
		
x.r.ns[j] <- mean(r.ns.total[2:6,j])	
				
		} #end of j

#####Hunting scenario
for(j in 4){
	n.s[1,1:8,j]<-N[y.end.sim,1:8]  #start with conditions at end of data
	n.s[1,9,j]<-0 
	n.s.total[1,j] <- N.total[y.end.sim]
	Re.s[1,j] <- Re[y.end.sim]
	F.s[1,j] <- F[y.end.sim]
	pv.cow.s[1,j] <-  pv.cow[y.end.sim]
	phi.s[1,j] <- phi[y.end.sim]
	calf.ratio.s[1,j] <- calf.ratio[y.end.sim] #calves
	cow.ratio.s[1,j] <- cow.ratio[y.end.sim]
	bull.ratio.s[1,j] <- bull.ratio[y.end.sim]
	delta.phi.s[1,j] <- 0
	r.ns[1,1,j]<-0
	r.ns[1,2,j]<-0  
	r.ns[1,3,j]<-0
	r.ns[1,4,j]<-0
	r.ns[1,5,j]<-0		
	r.ns[1,6,j]<-0
 	r.ns[1,7,j]<-0
	r.ns[1,8,j]<-0
	r.ns[1,9,j]<-0
	for(t in 2:y.scen.yrs+1){ #j indexes the scenario number
	
	#get random number of removals
	
	
	r.hunt[t]<-min(prop.remove[t]*n.s.total[t-1,j],1500)
	
	#removal vector
	r.ns[t,1,j]<-r.hunt[t] * n.s[t-1,1,j]/n.s.total[t-1,j]
	r.ns[t,2,j]<-r.hunt[t] * n.s[t-1,2,j]/n.s.total[t-1,j]  
	r.ns[t,3,j]<-r.hunt[t] * n.s[t-1,3,j]/n.s.total[t-1,j]
	r.ns[t,4,j]<-r.hunt[t] * n.s[t-1,4,j]/n.s.total[t-1,j]	
	r.ns[t,5,j]<-r.hunt[t] * n.s[t-1,5,j]/n.s.total[t-1,j] 
	r.ns[t,6,j]<-r.hunt[t] * n.s[t-1,6,j]/n.s.total[t-1,j]
	r.ns[t,7,j]<-r.hunt[t] * n.s[t-1,7,j]/n.s.total[t-1,j]
	r.ns[t,8,j]<-r.hunt[t] * n.s[t-1,8,j]/n.s.total[t-1,j]	
	r.ns[t,9,j]<-0 #class for vacination in later runs
	
	S[t,1,j] <- max(0,p[1]-p[1]^(1-q)*r.ns[t,1,j]/n.s[t-1,1,j])
	S[t,2,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,2,j]/n.s[t-1,2,j])
	S[t,3,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,3,j]/n.s[t-1,3,j])
	S[t,4,j] <- max(0,p[1]-p[1]^(1-q)*r.ns[t,4.,j]/n.s[t-1,4,j])
	S[t,5,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,5,j]/n.s[t-1,5,j])
	S[t,6,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,6,j]/n.s[t-1,6,j])
	S[t,7,j] <- max(0,p[2]-p[2]^(1-q)*r.ns[t,7,j]/n.s[t-1,7,j])
	S[t,8,j] <- max(0,p[3]-p[3]^(1-q)*r.ns[t,8,j]/n.s[t-1,8,j])
	S[t,9,j] <- 0
	
	
		
	#	
##frequency depdendent transmission
s_numer_phi[t,j]<-S[t,6,j]*n.s[t-1,6,j] + S[t,7,j]*psi*n.s[t-1,7,j]
s_denom_phi[t,j] <- ((1-m)*(S[t,1,j]*n.s[t-1,1,j]+S[t,4,j]*n.s[t-1,4,j]) + sum(S[t,2:3,j]*n.s[t-1,2:3,j]) + sum(S[t,5:7,j]*n.s[t-1,5:7,j]))^z
phi.s[t,j] <- 1-exp(-beta * s_numer_phi[t,j] /s_denom_phi[t,j])

#####Transition Matrix
#row 1
B[1,3,t,j]<-S[t,3,j]*f.neg[t]*(1-phi.s[t,j])
B[1,6,t,j]<-S[t,6,j]*f.con[t]*(1-v)*(1-phi.s[t,j])
B[1,7,t,j]<-S[t,7,j]* (f.pos[t] *(1-psi)*(1-phi.s[t,j])  +  f.pos[t] * psi*(1-v)*(1-phi.s[t,j]))
B[1,9,t,j]<-0


#row 2
B[2,1,t,j]<-S[t,1,j]*(1-m)*(1-phi.s[t,j])
#row 3
B[3,2,t,j]<-S[t,2,j]*(1-phi.s[t,j])
B[3,3,t,j]<-S[t,3,j]*(1-phi.s[t,j])

#row4
B[4,3,t,j]<-S[t,3,j]*f.neg[t]*phi.s[t,j]
B[4,6,t,j]<-S[t,6,j]*f.con[t]*(v+phi.s[t,j]-v*phi.s[t,j])
B[4,7,t,j]<-S[t,7,j]*(f.pos[t]*psi*(v+phi.s[t,j]-v*phi.s[t,j]) + f.pos[t]*(1-psi)*phi.s[t,j])
B[4,9,t,j]<-0
#row5
B[5,1,t,j]<-S[t,1,j]*phi.s[t,j]*(1-m)
B[5,4,t,j]<-S[t,4,j]*(1-m)

#row6
B[6,2,t,j]<-S[t,2,j]*phi.s[t,j]
B[6,3,t,j]<-S[t,3,j]*phi.s[t,j]
B[6,5,t,j]<-S[t,5,j]
B[6,7,t,j] <-0

#row 7
B[7,6,t,j]<-S[t,6,j]
B[7,7,t,j]<-S[t,7,j]

#row 8
B[8,1,t,j]<-S[t,1,j]*m
B[8,4,t,j]<-S[t,4,j]*m
B[8,8,t,j]<-S[t,8,j]


#row 9
B[9,9,t,j]<-S[t,9,j]
	
#estimates of mean of posterior distribution
mu.s[t,1:9,j] <- B[1:9,1:9,t,j]%*%n.s[t-1,1:9,j]

#get logs of deterministic predictions
for(k in 1:9){
	log_mu.s[t,k,j] <- log(max(1,mu.s[t,k,j]))
	}

log_n.s[t,1:9,j] ~ dmnorm(log_mu.s[t,1:9,j],tau.p*y.I2) 

#exponentiate log of population size
for(k in 1:9){
	n.s[t,k,j] <- min(exp(log_n.s[t,k,j]),10000)
	}

	
#derived quantities
Re.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / (n.s[t-1,6,j]+psi*n.s[t-1,7,j])
F.s[t,j] <- phi.s[t,j] * sum(S[t,1:3,j]*n.s[t,1:3,j]) / sum(n.s[t-1,1:3,j]) 
pv.cow.s[t,j] <-  sum(n.s[t,6:7,j])/(sum(n.s[t,6:7,j])+n.s[t,3,j])
n.s.total[t,j] <- sum(n.s[t,1:8,j])
cow.ratio.s[t,j] <- (n.s[t,3,j] + sum(n.s[t,6:7,j])) / n.s.total[t,j] #cows
calf.ratio.s[t,j] <- (n.s[t,1,j] + n.s[t,4,j]) / n.s.total[t,j] #calves
bull.ratio.s[t,j] <- n.s[t,8,j] / n.s.total[t,j] #bulls

delta.phi.s[t,j] <-phi.s[t,j] - phi[y.end.sim]


r.ns.total[t,j] <- sum(r.ns[t,1:8,j])

		} #end of t
		
x.r.ns[j] <- mean(r.ns.total[2:6,j])	
				
		} #end of j

	
####vacination##########################
##
for(j in 5){
	n.s[1,1:8,j]<-N[y.end.sim,1:8]  #start with conditions at end of data
	n.s[1,9,j]<-0 
	n.s[1,10,j]<-0
	n.s.total[1,j] <- N.total[y.end.sim]
	Re.s[1,j] <- Re[y.end.sim]
	F.s[1,j] <- F[y.end.sim]
	
	pv.cow.s[1,j] <-  pv.cow[y.end.sim]
	phi.s[1,j] <- phi[y.end.sim]
	calf.ratio.s[1,j] <- calf.ratio[y.end.sim] #calves
	cow.ratio.s[1,j] <- cow.ratio[y.end.sim]
	bull.ratio.s[1,j] <- bull.ratio[y.end.sim]
	delta.phi.s[1,j] <- 0
	
	for(t in 2:y.scen.yrs+1){ #j indexes the scenario number
	
	
	
	S[t,1,j] <- p[1]
	S[t,2,j] <- p[2]
	S[t,3,j] <- p[2]
	S[t,4,j] <- p[1]
	S[t,5,j] <- p[2]
	S[t,6,j] <- p[2]	
	S[t,7,j] <- p[2]
	S[t,8,j] <- p[3]
	S[t,9,j] <- p[2]
	S[t,10,j] <-p[2]
	
	
		
	phi.s[t,j] <- 1-exp(-beta *p[2]^q[1]* (n.s[t-1,6,j] + psi*n.s[t-1,7,j] + epsilon*n.s[t-1,10,j])   / ( (1-m)*p[1]^q[1]*(n.s[t-1,1,j] + n.s[t-1,4,j]) + p[2]^q[1]*(sum(n.s[t-1,2:3,j]) + sum(n.s[t-1,5:7,j]) + sum(n.s[t-1,9:10,j]) ))^z)

		
	
#####Transition Matrix
#row 1
B[1,3,t,j]<-S[t,3,j]*f.neg[t]*(1-phi.s[t,j])
B[1,6,t,j]<-S[t,6,j]*f.con[t]*(1-v)*(1-phi.s[t,j])
B[1,7,t,j]<-S[t,7,j]* (f.pos[t] *(1-psi)*(1-phi.s[t,j])  +  f.pos[t] * psi*(1-v)*(1-phi.s[t,j]))
B[1,9,t,j]<-S[t,9,j]*f.neg[t]*(1-phi.s[t,j])
B[1,10,t,j]<-S[t,10,j]*f.con[t]*(1-epsilon*v)*(1-phi.s[t,j])

#row 2
B[2,1,t,j]<-S[t,1,j]*(1-m)*(1-phi.s[t,j])
#row 3
B[3,2,t,j]<-S[t,2,j]*(1-phi.s[t,j])
B[3,3,t,j]<-S[t,3,j]*(1-phi.s[t,j])

#row4
B[4,3,t,j]<-S[t,3,j]*f.neg[t]*phi.s[t,j]
B[4,6,t,j]<-S[t,6,j]*f.con[t]*(v+phi.s[t,j]-v*phi.s[t,j])
B[4,7,t,j]<-S[t,7,j]*(f.pos[t]*psi*(v+phi.s[t,j]-v*phi.s[t,j]) + f.pos[t]*(1-psi)*phi.s[t,j])
B[4,9,t,j]<-0
B[4,10,t,j]<-S[t,10,j]*f.con[t]*(epsilon*v+phi.s[t,j]-epsilon*v*phi.s[t,j])


#row5
B[5,1,t,j]<-S[t,1,j]*phi.s[t,j]*(1-m)
B[5,4,t,j]<-S[t,4,j]*(1-m)

#row6
B[6,2,t,j]<-S[t,2,j]*phi.s[t,j]
B[6,3,t,j]<-S[t,3,j]*phi.s[t,j]
B[6,5,t,j]<-S[t,5,j]
B[6,7,t,j] <-0

#row 7
B[7,6,t,j]<-S[t,6,j]
B[7,7,t,j]<-S[t,7,j]
B[7,10,t,j] <- p[2]

#row 8
B[8,1,t,j]<-S[t,1,j]*m
B[8,4,t,j]<-S[t,4,j]*m
B[8,8,t,j]<-S[t,8,j]


#row 9
B[9,9,t,j]<-S[t,9,j]*(1-phi.s[t,j])


#row 10
B[10,9,t,j] <-p[2]*phi.s[t,j]
r.remove[t,j] <-min(prop.remove[t]*n.s.total[t-1,j],1500)

#vector for removals
	#removal vector
	r.ns[t,1,j]<-0
	r.ns[t,2,j]<-r.remove[t,j]*n.s[t-1,2,j]/n.s.total[t-1,j]	
	r.ns[t,3,j]<-r.remove[t,j]*n.s[t-1,3,j]/n.s.total[t-1,j] 
	r.ns[t,4,j]<-0
	#seropositives removed
	r.ns[t,5,j]<-0
	r.ns[t,6,j]<-0
	r.ns[t,7,j]<-0
	r.ns[t,8,j]<-0
	r.ns[t,9,j]<- -(r.ns[t,2,j] + r.ns[t,3,j]) #This quantity is negative so that it will be added in the model for mu below 
	r.ns[t,10,j]<-0

	
#estimates of mean of posterior distribution
mu.s[t,1:10,j] <- B[1:10,1:10,t,j]%*%n.s[t-1,1:10,j] - r.ns[t,1:10,j]



#get logs of deterministic predictions
for(k in 1:10){
	log_mu.s[t,k,j] <- log(max(1,mu.s[t,k,j]))
	}

log_n.s[t,1:10,j] ~ dmnorm(log_mu.s[t,1:10,j],tau.p*y.I3) 

#exponentiate log of population size
for(k in 1:10){
	n.s[t,k,j] <- min(exp(log_n.s[t,k,j]),10000)
	}
	

#derived quantities
Re.s[t,j] <- (phi.s[t,j] * p[2] * sum(n.s[t-1,2:3,j]) + phi.s[t,j] * p[2] * n.s[t-1,10,j] + phi.s[t,j] * p[1] * n.s[1,t-1,j] + n.s[t,4,j] ) / (n.s[t-1,6,j] + n.s[t-1,10,j] + psi*n.s[t-1,7,j])

F.s[t,j] <- (phi.s[t,j] * p[2] * sum(n.s[t-1,2:3,j]) + phi.s[t,j] * p[2] * n.s[t-1,10,j] + phi.s[t,j] * p[1] * n.s[1,t-1,j] + n.s[t,4,j] )  / (sum(n.s[t-1,1:3,j])+n.s[t-1,9,j])

pv.cow.s[t,j] <-  (sum(n.s[t,6:7,j]) + n.s[t,10,j]) / (sum(n.s[t,6:7,j]) + n.s[t,3,j] + sum( n.s[t,9:10,j] ))

n.s.total[t,j] <- sum(n.s[t,1:10,j])
cow.ratio.s[t,j] <- (n.s[t,3,j] + sum(n.s[t,6:7,j] + sum(n.s[t,9:10,j]))) / n.s.total[t,j]
calf.ratio.s[t,j] <- (n.s[t,1,j] + n.s[t,4,j]) / n.s.total[t,j] #calves
bull.ratio.s[t,j] <- n.s[t,8,j] / n.s.total[t,j] #bulls
delta.phi.s[t,j] <-phi.s[t,j] - phi[y.end.sim]
r.ns.total[t,j] <- sum(r.ns[t,1:8,j])

		} #end of t
		
x.r.ns[j] <- mean(r.ns.total[2:6,j])	
				
		} #end of j
}
/*End of Model*/